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_L Abstract 

We develop an optimization-based approach to the problem of reconstructing 
c/5 temperature-dependent material properties in complex thermo-fluid systems described 

by the equations for the conservation of mass, momentum and energy. Our goal is to 
^ estimate the temperature dependence of the viscosity coefficient in the momentum 

equation based on some noisy temperature measurements, where the temperature is 
, governed by a separate energy equation. We show that an elegant and computation- 

ally efficient solution of this inverse problem is obtained by formulating it as a PDE- 
constrained optimization problem which can be solved with a gradient-based descent 
method. A key element of the proposed approach, the cost functional gradients are 
CN characterized by mathematical structure quite different than in typical problems of 

^ PDE-constrained optimization and are expressed in terms of integrals defined over the 

level sets of the temperature field. Advanced techniques of integration on manifolds 
are required to evaluate numerically such gradients, and we systematically compare 
three different methods. As a model system we consider a two-dimensional unsteady 
^ flow in a lid-driven cavity with heat transfer, and present a number of computational 

t> tests to validate our approach and illustrate its performance. 

>< 

^ Keywords: parameter estimation, material properties, optimization, adjoint anal- 

ysis, integration on level sets. 

1 Introduction 

In this work we propose and validate a computational approach to the reconstruction of 
material properties in complex multiphysics phenomena based on incomplete and possibly 
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noisy measurements. The material properties we are interested in here are the transport co- 
efficients characterizing diffusion processes such as the viscosity or the thermal conductivity, 
and we focus on problems in which these coefficients depend on the state variables in the 
system. By the "multiphysics" aspect we mean situations in which the material property 
used in one conservation equation is a function of a state variable governed by a different 
conservation equation, e.g., reconstruction of the temperature dependence of the viscosity 
coefficient used in the momentum equation, where the temperature is governed by a sepa- 
rate energy equation, which is the specific model problem investigated in this study. This 
research is motivated by questions arising in the computational analysis and optimization of 
advanced welding processes which involves modelling complex alloys in the liquid phase at 
high temperatures [l]. 

Inverse problems of parameter estimation for partial differential equations (PDEs) have 
received significant attention in the literature, both regarding theoretical |2| and practical 
aspects [3j. However, most of the attention focused on problems in which the material prop- 
erties are functions of the space variable (i.e., the independent variable in the problem). 
Such problems are, at least in principle, relatively well understood and represent the foun- 
dation of, for example, numerous imaging techniques in medicine ^ and earth sciences 
The problem considered here is in fact different, in that the material properties are sought 
as functions of the state (dependent) variables in the system which gives rise to a number 
of computational challenges absent in the "classical" parameter estimation problem. Other 
than the seminal work of Chavent and Lemonnier [6j, our earlier study [7| concerning a 
simplified model problem and a few investigations of fully discrete formulations surveyed 
in [7], there does not seem to be much literature concerning computational methods for 
this type of parameter estimation problems. One way to solve such inverse problems is to 
formulate them as optimization problems and this is the approach we will follow focusing on 
the "optimize-then-discretize" paradigm in which the optimality conditions are formulated 
at the continuous (PDE) level and only then discretized. The goal of this investigation is to 
extend the approach formulated in [T] for a simple model to a realistic multiphysics problem 
involving time-dependent fluid flow in a two-dimensional (2D) domain. As will be shown 
below, a number of computational difficulties will need to be overcome in order to achieve 
this goal. 

As a key contribution of this work, we address a number of computational challenges 
related to accurate and efficient evaluation of cost functional gradients which are critical to 
the implementation of the proposed approach. More specifically, these gradients are given in 
terms of integrals of expressions involving state and adjoint variables defined on a grid over 
contours given by the level sets of the temperature field. A number of techniques have been 
proposed for the numerical evaluation of integrals defined over manifolds defined by level-set 
functions. Some of them rely on regularized Dirac delta and Heaviside functions IsIIq], or 



discretization of the Dirac delta function 10 12 . Similar approaches, based on approxima- 



tions of the Dirac delta functions obtained using the level-set function and its gradient, were 



developed by Towers 13,14 . The family of geometric approaches proposed by Min and Gi- 



bou in 15 , 16 relies on a decomposition of the entire domain into simplices. We emphasize 
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that the problem discussed in this work is in fact more complicated, as the computation 
of our cost functional gradients requires evaluation of the corresponding integrals for the 
level-set values spanning the entire state space of interest, hence there are also additional 
issues related to the discretization of the state space which were outside the scope of refer- 
ences |8-16 . Thus, in order to address these questions and assess the different trade-offs in 



the choice of the numerical parameters we will compare the computational performance of 
three different methods for the evaluation of cost functional gradients. 

The structure of this paper is as follows: in the next Section we formulate our model 
problem, in the following Section we cast the problem of parameter estimation as an opti- 
mization problem, an adjoint-based gradient-descent algorithm is formulated in Section |4| 
in Section [5] we outline some regularization strategies needed in the presence of measurement 
noise, whereas in Section [6] we analyze three different numerical approaches to the evaluation 
of the cost functional gradients; extensive computational results are presented in Section [7] 
with discussion and conclusions deferred to Section M 



2 Model Problem 

Let Q C M"*, d = 2,3, be the spatial domain on which our model problem is formulated. 
To fix attention, but without loss of generality, in the present investigation we focus on 
the problem of a reconstruction of the temperature dependence fi = yu(T) of the viscosity 
coefficient : M — )■ M"*" in the momentum equation (Navier-Stokes equation), where the 
temperature T is governed by a separate energy equation 

dtu + (u • V)u + Vp - V ■ [/i(T) [Vu + ( Vu)^]] =0 in Q, (la) 

V ■ u =0 inQ, (lb) 
dtT + (u • V)T - V • [kVT] =0 in fi, (Ic) 

subject to appropriate Dirichlet (or Neumann) boundary and initial conditions 

u = Us on (9fi, (2a) 

T = Tb on dn, (2b) 

u(-,0) =uo, T(-,0) =ro infi. (2c) 

The specific inverse problem we address in this investigation is formulated as follows. Given 
a set of time-dependent "measurements" {Ti(t)}fii of the state variable (temperature) T 
at a number of points {^i}fli in the domain Q (or along the boundary dQ) and obtained 
within the time window t G [0,t/], we seek to reconstruct the constitutive relation fi = fi{T) 
such that solutions of problem ([l])-([2]) obtained with this reconstructed function will fit best 
the available measurements. 

In regard to reconstruction of constitutive relations in general, it is important that such 



relations be consistent with the second principle of thermodynamics 17 . There exist two 



mathematical formalisms, one due to Coleman and Noll 18 and another one due to Liu 119 
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developed to ensure in a very general setting that a given form of the constitutive relation 
does not violate the second principle of thermodynamics. In continuous thermodynamical 
and mechanical systems this principle is expressed in terms of the Clausius-Duhem inequality 
[20] which in the case of our present model problem ([T])-([2]) reduces to the statement that 
/x(T) > for all possible values of T. 

In our discussion below we will also need definitions of the following intervals, cf. Figure 

• [Tq,, Tp] = [min^gf^ T{x), max^g^^ ^(x)] which represents the temperature range spanned 
by the solution of problem ([T]); thus, following 21 , we will refer to the interval 
X = [Tq,T^] as the identifiability interval, 

• L = [Ta,Tf,], where < and T^, > Tp] this will be the temperature interval on 
which we will seek to obtain a reconstruction of the material property; we note that in 
general the interval L will be larger than the identifiability interval, i.e., X C L, and 

• M. = [mini<j<M mino<t<f^ Tj(t), maxi<j<A/ maxo<t<t^ Tj(t)] which defines the temper- 
ature range spanned by the measurements {Ti}fi^^\ this interval is always contained 
the identifiability interval X, i.e., C X; we will refer to the interval Ai as the 
measurement span. 



3 Parameter Estimation as an Optimization Problem 

It is assumed that the constitutive relations /^(T) are differentiable functions of the state 
variable (temperature) and belong to the following set 

= {/i(T) piecewise on L; < < fi{T) < < 00, VT G L}, (3) 

where m^, G M"*". We will also assume that the set consisting of constitutive relations 
/i(T) defined on L is embedded in a Hilbert (function) space X to be specified below. Solving 
our parameter estimation problem is therefore equivalent to finding a solution to the operator 
equation 

J-(/i) = T, (4) 

where J-" : — )■ (L2([0, t/]))^ is the map from the constitutive relations to the measure- 
ments. An approach commonly used to solve such problems consists in reformulating them 
as least-squares minimization problems which in the present case can be done by defining 
the cost functional J : — )■ M as 

M 
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Figure 1: Schematic showing (left) the solution T(to, x) at some fixed time to and (right) 
the corresponding constitutive relation /i(T) defined over their respective domains, i.e., Q = 
(— 1, 1) and the identifiability region X. The thick dotted line represents an extension of the 
constitutive relation yu(T) from X to the interval L. In the Figure on the right the horizontal 
axis is to be interpreted as the ordinate. 
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where the dependence of the temperature field T{-; fi) on the form of the constitutive relation 
/i = /i(T) is given by governing system ([T])-([2]). The optimal reconstruction p, is obtained as 
a minimizer of cost functional ([s]), i.e., 

/t = argmin j7'(/u). (6) 

We recall that the constitutive property is required to satisfy the positivity condition /i(T) > 
for all T G L, cf. ([s]). Therefore, the optimal reconstruction fi should in fact be obtained 
as an inequality-constrained minimizer of cost functional ([s]), i.e., 

/t = argmin (7) 



(T)>o, TeL 



We add that in problems involving constitutive relations depending on several state vari- 
ables the inequality constraint > will be replaced with a more general form of 
the Clausius-Duhem inequality ^O]. Different computational approaches for converting 
inequality-constrained optimization problems to unconstrained formulations are surveyed 
in (22 ,23 . Here we follow a straightforward approach based on the so-called "slack" vari- 



able [24]. We define a new function 6{T) : M -> M such that 

fi{T)=e\T)+m^, 



where is a lower bound for jJ.{T), cf. (|3]). This change of variables allows us to transform 
the inequality-constiained optimization problem ([T]) to a new unconstrained one 

^ = argmin J^(^), (9) 

where the constraint /i(T) > is satisfied automatically when minimization is performed 
with respect to the new variable 0{T). In view of ([s]), we note that the new optimization 
variable 9 belongs to the following set 



Se = {9{T) piecewise on L; \e{T)\ < ^M^ - m^, VT G L}. (10) 
The governing PDE system ^ can thus be rewritten in the form 
dtM + (u • V)u + Vp - V ■ W{T) + m^)[Vu + (Vu" " 



Vu 



redefining cost functional ([s]) in terms of the new variable 

1 f'f 



=0 


in fi. 


(11a) 


=0 


in fi. 


(lib) 


=0 


in fi. 


(11c) 


new 


problem ^ 


requires 



i=i 
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Problem (|9| is characterized by the first-order optimahty condition which requires the 
Gateaux differential of cost functional (12), defined as J'{6]6') = hra^^oe^^[J'(9 + e9') — 



J'{0)], to vanish for all perturbations 6' G X |25j, i.e., 

ye'ex j'{e;e') = o. 



(13) 



The (local) optimizer 6 can be computed with the following gradient descent algorithm as 
9 = lim„_^oo where 



n = 1, 



(14) 



in which ^eJ'{9) represents the gradient of cost functional x7{9) with respect to the control 
variable 9 (we will adopt the convention that a subscript on the operator V will be used when 
differentiation is performed with respect to variables other than x), the length of the 

step along the descent direction at the n-th iteration, whereas = y/fJ'O — "'^^ is the initial 



guess taken, for instance, corresponding to a constant /io, or some other appropriate initial 
approximation. For the sake of clarity, formulation (14) represents the steepest-descent 



algorithm, however, in practice one typically uses more advanced minimization techniques, 
such as the conjugate gradient method, or one of the quasi-Newton techniques 26 . We note 



that, since minimization problem ^ is in general nonconvex, condition (13) characterizes 
only a local, rather than global, minimizer. 



The key ingredient of minimization algorithm ( 14 ) is computation of the cost functional 
gradient '^eJ{G)- We emphasize that, since 9 = 9{T) is a continuous variable, the gradient 
'^eJ{9) represents in fact the infinite-dimensional sensitivity of J{9) to perturbations of 
9{T). This gradient can be determined based on suitably defined adjoint variables (Lagrange 
multipliers) obtained from the solution of the corresponding adjoint system. Since this 
derivation differs in a number of imported technical details from analogous derivations in 
"standard" PDE-constrained optimization problems, it will be reviewed in Section |4j The 



expression for the gradient is then validated for consistency in Section 7.4 



4 Cost Functional Gradients via Adjoint— based Anal- 
ysis 



Since the new variable 9{T) belongs to set Se, cf. (10), we will seek to reconstruct 9{T) as 
elements of the Sobolev space H^{L), so that the gradient V^JT will need to be obtained 
with respect to the corresponding inner product. However, in order to make the derivation 
procedure easier to follow, we will first obtain an expression for the gradient in the space 
L2(L), and only then will obtain the Sobolev gradients which will be eventually used in the 
solution of optimization problem ([9]). In all these steps our transformations will be formal 



We begin by computing the directional (Gateaux) differential of cost functional (12) which 
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yields 



J'ie-e') 



fi{r)]T'{T,^i-e,e') dr, 



(15) 



where the perturbation variable T'{6, 9') satisfies the perturbation system obtained from 
([T])-(|2]). Next, we invoke the Riesz representation theorem 27 for the directional differential 
J\e; ■) which yields 



J'ie-d') 



X 



(16) 



where (■, ■)x represents an inner product in the Hilbert space X (we will first set X = L2(L) 
and afterwards change this to X = H^{L)). We note that the expression on the right-hand 
side (RHS) in (15) is not consistent with Riesz representation (16), since, as will be shown 
below, the perturbation variable 6' is hidden in the system defining T'{6,6'). However, this 
expression can be transformed to Riesz form (16) with the help of a suitably-defined adjoint 
variable, a result which is stated in Theorem 4.1 below. The main aspect in which this 
derivation differs from the standard adjoint analysis |28] is that the inner product in Riesz 
identity (16) is defined using the state variable (temperature) as the integration variable, 
whereas the variational formulation of ([T])-([2]) is defined using integration with respect to 
the independent variables (x and t). 

Theorem 4.1. Let Q be a sufficiently regular open hounded domain and 9' & X = H^{L). 
We assume that the solutions u and T of system ([T]) -([2]) are sufficiently smooth. Then, the 
Riesz representation of directional differential (|T5|) has the form 



J'{e;e') 



5(T(x)-.)0(.) 



oo J fl 



[Vu + (Vu)^] : Vu* dr 



^'(s)dxds, (17) 



where 6{-) denotes Dirac delta function and the adjoint state {u*,T*} is defined as the 
solution of the system 



-dtu* 



(u ■ V)u* - V ■ 0-* + u* ■ (Vu) 



T*vr = 



in f2, 



V ■ u* 



de 



-dtT* - (u ■ V)T* - V ■ [kVT*] + 29{T) ^(T) [Vu + (Vu)"^] : V*u 



A/ 



^[T(x,;0)-f,]5(x-x, 



1=1 



:i8a) 



in Q, 



;i8b) 



in 



:i8c) 



where a* = -p*X+[e'^(T)+m^) [V u* + (Vu*)"^] , with the following boundary and terminal 
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conditions 



u 






0, T*(-:t 



on dQ, 
on dQ, 
in Q. 



(19) 



Proof. We will denote the stress tensor cr = —pX+ {6'^{T)+m^) [Vu + (Vu)-^]] and rewrite 
the governing system (11) as 



dtVL + (u ■ V)u - V ■ cr = in fi, 
V ■ u = inVL, 
dtT + (u • V)T - V ■ [fcVT] =0 in n 



(20) 



with the boundary and initial conditions (j2]). Perturbing the state variables u, p and T, 
which are functions of time and space, we get 



u = uo + eu' + 0(e^), 
p = Po + ep' + 0(e^), 
T = To + eT' + 0(e2), 



(21) 



so that the corresponding expansion of the constitutive relation 6{T) will have the following 
form 



9{T) = ^o(T) + ee'iT) + 0(e^) = ^o(To) + e^(To) r(0o; ^^') + e^'(To) + 0(e 



(221 



where the subscript "0" is used to denote the unperturbed (reference) material property and 
the state variable, whereas the prime denotes the corresponding perturbations. We also have 



de 



d\T) = ^^o'(To) + 2e9,{To) —{To) T'iOo; 9') + 2e^^o(To)e'(To) + &{e' 



(23) 



Substituting (21) and (23) into (20), collecting terms corresponding to e in different powers 
and denoting 



~ A 

cr = 



{2e{T) ^(T) T\9o; e') + 29{T) 9'{T)^ [Vu + (Vu)^ 



we obtain the perturbation (sensitivity) system corresponding to (11) 



(9tu' + (u' • V)u + (u • V)u' - V • (6- + a-) =0 in fi, 

V ■ u' =0 in 
dtT' + (u ■ V)T + (u ■ V)T' - V • [kVT'] =0 in Vt 



(24a) 
(24b) 
(24c) 
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with the following boundary and initial conditions 



u' = on dfl, 

T' = on dn, 

u (-,0) = 0, r'(-,0) = inn. 



(25a) 
(25b) 
(25c) 



Then, integrating equation (24a) against u*, equation (24b) against p*, and equation ( |24c ) 
against T* over the space domain Q and time [0,t/], integrating by parts and factorizing u', 
T' and p', we arrive at the following relation 



[-dtu* + u* ■ ( Vu)^ - (u ■ V)u* - V ■ (7* + T* VT] ■ u' dr 
(V • u* )p'dxdT 



Jn 



Jn 



+ 
+ 



Jn I 



dB 

-dtT* - (u ■ V)T* - V ■ (kVT*) + 2e{T) — {T)[Vu + {Vuf] : Vu* 

(Jj-L 



T' d:si dr 



Jn 



2^(T)0'(T)[Vu+ (Vu)^] : Vu*dxdr = 0. 



(26) 

We now require that the adjoint variables u*, p* and T* satisfy system (18)-(19). We also 
note that owing to the judicious choice of the RHS term in (18c), the last term in relation 
(26) is in fact equal to the directional differential J''{6;6'), so that we have 

'0 Jn 



0(T(x, r))0'(T(x, r)) [Vu(x, r) + (Vu(x, T)f] : Vu*(x, r) rfxc/r. 



(27) 

where, for emphasis, we indicated the integration variables as arguments of the state and 
adjoint variables. We note that this expression is still not in Riesz form (16), where inte- 
gration must be performed with respect to the state variable (temperature T). Thus, we 
proceed to express for any given function /(T) its pointwise evaluation at T(x) through the 
following integral transform. It is defined using a "change-of-variable" operator, denoted 11, 
such that for given functions / : M — )■ M and T : — M, we have 



/+00 
5(T(x)-.)/(.)t/.4(n/)(x). 
-oo 



(28) 



Using this transform to express /(T(x)) = 0(T(x, r))0'(T(x, r)) in (27) and changing the 
order of integration (Fubini's Theorem), we obtain expression (17) which is the required 
Riesz representation (16) of directional differential (15). □ 

We remark that we were able to prove an analogous result using a simpler approach based 
on the Kirchhoff transform in [T] , where both the constitutive relation and the state variable 
were governed by the same equation (i.e., the problem was not of the "multiphysics" type). 
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With the Riesz representation estabhshed in (17), we now proceed to identify expressions 
for the cost functional gradient qJ according to (16) using different spaces X. While this 
is not the gradient that we will use in actual computations, we analyze first the "simplest" 
case when X = L2(L), i.e., the space of functions square integrable on [Ta,Tb], as it already 
offers some interesting insights into the structure of the problem. The L2 gradient of the 
cost functional hence takes the form 



'^^J{s) = -2 r [ 6{T{ 
Jo Jn 



[Vu + (Vu 



Vu* (hi-dr. 



(29) 



As was discussed at length in |7|, the L2 gradients are not suitable for the reconstruction of 
material properties in the present problem, because in addition to lacking necessary smooth- 
ness, they are not defined outside the identifiability region (other than perhaps through 
a trivial extension with zero). Given the regularity required of the constitutive relations, 
(10), the cost functional gradients should be elements of the Sobolev space -ff^(L) of 

i/i(L), we 



cf 

functions with square-integrable derivatives on L 
obtain 



Using (16), now with X 



J'{e;e') 



J, e' 



Ta 



ds ds 



(30) 



ds 



[we note that the L2 
Performing integration by parts with the 



in which £ G M is a parameter with the meaning of a "temperature-scale'' 
inner product is recovered by setting £ = in (30) 
assumption that the Sobolev gradient J7 satisfies the homogeneous Neumann boundary 
conditions at T = Ta,T;,, and noting that relation (30) must be satisfied for any arbitrary 
6', we conclude that the Sobolev gradient can be determined as a solution of the following 
inhomogeneous elliptic boundary- value problem where the state variable (temperature) acts 
as the independent variable 



^2 

«2 " -nH' 



ds 

d 

ds 



Vfj = 



for s = Ta,Tb. 



(31a) 
(31b) 



We recall that by changing the value of the temperature-scale parameter i we can control 
the smoothness of the gradient i7(^); and therefore also the relative smoothness of the 
resulting reconstruction of 0{T), and hence also the regularity of yu(T). More specifically, as 
was shown in 29 , extracting cost functional gradients in the Sobolev spaces H^, p > 0, is 
equivalent to applying a low-pass filter to the L2 gradient with the quantity i representing 
the "cut-off" scale. There are also other ways of defining the Sobolev gradients in the 
present problem which result in gradients characterized by a different behavior outside the 
identifiability region. These approaches were thoroughly investigated in [t], and since they 
typically lead to inferior results, they will not be considered here. We finally conclude that 
iterative reconstruction of the constitutive relation /i(T) involves the following computations 



11 



1. solution of direct problem (11) with boundary and initial conditions ([2]), 



2. solution of adjoint problem (18)-(19) 



3. evaluation of expression (29) for the cost functional gradient, 



4. computation of the smoothed Sobolev gradient by solving (31) 



While steps (1), (2) and (4) are fairly straightforward, step (3) is not and will be thoroughly 
investigated in Section [6] 

As we discussed in detail in [t], while the Sobolev gradient J may be defined on 
an arbitrary interval L D X, the actual sensitivity information is essentially available only 
on the identifiability interval X (see Figure [T]). In other words, extension of the gradient 



outside X via (31) does not generate new sensitivity information. Since, as demonstrated 
by our computational results reported in [T], such techniques are not capable of accurately 
reconstructing the relation /i(T) on an interval L much larger than the identifiability region 
X, here we mention a different method to "extend" the identifiability region, so that the re- 
lation yu(T) can be reconstructed on a larger interval. This can be done in a straightforward 
manner by choosing appropriate time-dependent boundary conditions for temperature Tb 



in (2b) which will result in a suitable identifiability region X and the measurement span M.. 
Computational results illustrating the performance of our approach with different identifia- 
bility regions obtained using this method will be presented in Section 7.5[ We remark that 



extending the identifiability region in this way is not possible in time-independent problems 
where an iterative approach has to be used involving solution of a sequence of reconstruction 
problems on shifted identifiability regions [7|. 

5 Reconstruction in the Presence of Measurement Noise 

In this Section we discuss the important issue of reconstruction in the presence of noise in the 
measurements. As can be expected based on the general properties of parameter estimation 
problems and as will be confirmed in Section 7.6, incorporation of random noise into the 



measurements leads to an instability in the form of small-scale oscillations appearing in the 
reconstructed constitutive relations. In the optimization framework a standard approach to 



mitigate this problem is Tikhonov regularization 30 in which original cost functional (12) 



is replaced with a regularized expression of the form 

Jx{Q)=J{Q)^\\Q-Qty^^y (32) 

where A G is an adjustable regularization parameter, ^(T) represents a constitutive 
relation which our reconstruction ^(T) should not differ too much from, whereas || ■ \\y{T) is the 
Hilbert space norm in which we measure the deviation {6 — 6). Thus, the regularization term 



in (32), i.e., the second one on the RHS, involves some additional information which needs 



to be specified a priori, namely, the choice of the reference relation 6{T) and the space 3^(X). 
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As regards the reference function 0{T), one natural possibility is to consider a constant value 
corresponding to a constant material property, and this is the solution we will adopt below. 
We recall here that 0{T) is in fact a "slack" variable and is related to the actual constitutive 
relation via ([s]). As regards the choice of the space ^(X), we will follow the discussion in 
and consider a regularization term involving derivatives, namely 3^(X) = H^{X), where H^{I) 
denotes the Sobolev space equipped with the semi-norm 



the regularization term in (32) then becomes 



\m{i) 



A 
2 



9-0 



fde 



ds 



rp \ds ds 

yielding the following L2 gradient of the regularized cost functional 



(S) ds^ V.,HKX); 



(33) 



dd 



5(r(x) - s) e{s) [Vu + (Vu)"^] : Vu* d^dr 

d^e' 



+ A <! - [6is - Tp) - 5{s - T,)] 



(34) 



We remark that in obtaining (34) integration by parts was applied to the directional deriva- 



tive of the regularization term. Expression (|34|) can now be used to compute the Sobolev 

ti- 



gradients as discussed in Section |4| We add that penalty term (33) is defined on the identi- 
fiability interval X which is contained in the interval L on which the Sobolev gradients are 
computed. Computational tests illustrating the performance of the Tikhonov regularization 
on a problem with noisy data will be presented in Section |7.6[ In that Section we will also 



briefly analyze the effect of the regularization parameter A. We add that the stability and 
convergence of Tikhonov regularization using the Sobolev norm in the regularization 
term and applied to an inverse problem with similar mathematical structure, but formulated 
for a simpler PDE than ([T]), was established rigorously in 21 . 



6 Numerical Approaches to Gradient Evaluation 

Without loss of generality, hereafter we will focus our discussion on the 2D case. For some 
technical reasons we will also assume that 



Vte[o,t,] meas{xGfi, |VT(t,x)| = 0} = 0, 



(35) 



i.e., that the temperature gradient may not vanish on subregions with finite area. (This 
assumption is naturally satisfied when the temperature evolution is governed by an equation 



of the parabohc type such as (Ic) 



A key element of reconstruction algorithm (14) is evaluation of the cost functional gra- 



dients given, in the L2 case, by expression (29). The difficulty consists in the fact that at 



every instant of time t G [0,t/] and for every value of s (i.e., the dependent variable), one 
has to compute a line integral defined on the level set 



^,(t)^{xG^], T(t,x) = s} 



(36) 
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of the temperature field T(t,x). To focus attention on the main issue, the time dependence 
will be omitted in the discussion below. The integrand expression in such integrals is given in 
terms of solutions of the direct and adjoint problems ([T|-(|2]) and (18)-(19) which are approx- 
imated on a grid. As will be shown below, this problem is closely related to approximation 
of one-dimensional Dirac measures in M'', an issue which has received some attention in the 
literature [8 16 . We will compare different computational approaches to this problem, and 
in order to better assess their accuracy, we will first test them on the generic expression 



fis) 



^(0(s,x))5((x) rfx 



(37) 



for which the actual formula for the cost functional gradient (29) is a special case (except 
for the time integration). In (37) the function 0(s,x) : M x — )■ M represents the field 
whose s- level sets define the contours of integration Tg, whereas the function g{x.) : i7 — t- M 
represents the actual integrand expression. We note that by setting 0(s,x) = T(x) — s, 
g{x) = 26'(T(x))[Vu(x) + (Vu(x))"^] : Vu*(x) and adding time integration in (37), we 



recover the original expression (29) for the cost functional gradient. We emphasize, however. 



that the advantage of using (37) with some simple, closed-form expressions for i;A(s,x) and 
g{x.) as a testbed is that this will make our assessment of the accuracy of the proposed 
methods independent of the accuracy involved in the numerical solution of the governing 
and adjoint PDEs (needed to approximate u, T and u*). 



In anticipation of one of the proposed numerical approaches, it is useful to rewrite (37) 
explicitly as a line integral 



da 



(38) 



which is valid provided |V0| for every x G Tg, cf. assumption (35) (proof of the 
equivalence of expressions (37) and (38) may be found, for example, in |31|). Formula 



(38) makes it clear that for a fixed value of s expression (29) for the cost functional gradient 



can be interpreted as a sum of line integrals defined on the instantaneous s-level sets of the 
temperature field T(t,x). 



The problem of accurate numerical evaluation of the expressions given by either (37) or 



(38) has received much attention, especially since the invention of the level-set approach 
by Osher and Sethian 32 . Traditionally, the problem of integration over codimension-1 
manifolds defined by a level-set function 0(x) is studied in terms of the numerical evaluation 
of either the left-hand side (LHS) or right-hand side (RHS) expression in the following 



relation, analogous to (37)-(38), 



(x)=0 



h{-x)da= / 5(0(x))|V0(x)|/i(x)rfx, 



(39) 



where, by absorbing the factor | V0(x)|~ into the definition of the function h : — )■ M, one 
bypasses the problem of the points x G Fs where |V0(x)| 
fall into two main groups: 



0, cf. (35). These approaches 
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(A) reduction to a line (contour) integral, cf. (38 



), or the LHS of (39), and 



(B) evaluation of an area integral, cf. (37), or the RHS of (39). 
In the context of this classification, the methods of geometric integration developed by Min 



and Gibou 15 , 16 fall into the first category. This approach is based on decomposing the 
domain Q into simplices, which in the simplest 2D case can be achieved via a standard tri- 
angulation, and then approximating the level sets given by 0(x) = with piecewise splines 



inside each simplex. Expression ( 38 ) then breaks up into a number of definite integrals which 
can be evaluated using standard quadratures. 

In practice, however, area integration techniques (B) seem to have become more popular. 
One family of such techniques relies on regularization 6^ of the Dirac delta function with a 
suitable choice of the regularization parameter e which characterizes the size of the support. 
While in the simplest case in which the parameter e is determined based on the mesh size the 
error is 0(1) l8j, recently developed approaches l8l|9] achieve better accuracy by adjusting e 



based on the local gradient | V0| of the level-set function. Another family of area integration 



approaches is represented by the work of Mayo 10 further developed by Smereka 11 where 
a discrete approximation 6 of the Dirac delta function was obtained. This approach can also 
be regarded as yet another way to regularize delta function 5(0(x)) using a fixed compact 
support in the one-dimensional (ID) space of values 0(x). In the second group of approaches 



we also mention consistent approximations to delta function obtained by Towers in 13 , 14 
using the level-set function and its gradient computed via finite differences. 



In our present reconstruction problem, we have to evaluate the gradient expression (29) 
for the whole range of T G L, hence the discretization of the interval L will also affect the 
overall accuracy of the reconstruction, in addition to the accuracy characterizing evaluation 
of the gradient for a particular value of T. This is an aspect of the present problem which is 
outside the scope of earlier investigations concerning evaluation of the contour integrals of 
grid-based data [8 - 16 . Thus, we need to understand how the interplay of the discretizations 
of the physical space fl with the step size h and the state space L with the step size Ht 
affects the accuracy of the reconstruction. In principle, one could also consider the effect 
of discretizing the time interval [0,t/], however, the corresponding step size is linked to h 
via the CFL condition, hence this effect will not be separately analyzed here. There are 
also questions concerning the computational complexity of the different approaches. We will 



consider below the following three methods to evaluate expression (37), or equivalently (38) 



which are representative of the different approaches mentioned above 
#1 line integration over approximate level sets which is a method from group A based on a 



simplified version of the geometric integration developed by Min and Gibou in 15 , 16 , 



#2 approximation of Dirac delta measures developed by Smereka in 11 which is an exam- 
ple of a regularization technique and utilizes the area integration strategy from group 
B, and 

#3 approximation of contour integrals with area integrals, a method which also belongs 
to group B and combines some properties of regularization and discretization of Dirac 
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delta measures discussed above l8l|9,ll ; more details about this approach, including 



an analysis for its accuracy, are provided in Section 6.3 



To fix attention, we now introduce two different finite-element (FEM) discretizations of 
the domain Q based on 



triangular elements fl^ , i = 1, . . . , N^, such that 



fi = |Jnf, and 



(40) 



i=l 



quadrilateral elements flf, i = 1, . . . , N^, such that 



□ 

i 1 



(41) 



i=l 



where N^^ and A^n are the total numbers of the elements for each type of discretization (in 
case of uniform triangulation one has A^a = 2A^n)- In our computational tests we will assume 
that the functions 0(s,x) and g{x) are given either analytically, or in terms of the following 
FEM representations 



■5,X) OA 



z = l,...,iVA, (42) 

= Sti^^^Ux), 9i^)\n^ = SLi^?fc^Ux), z = 1, . . . ,iVa, (43) 
where 0^ and are the given nodal values of the functions 0(s, x) and 5'(x), whereas ipK'x.) 



are the basis functions (linear in (42) and bilinear in (43) [33]). We also discretize the 
reconstruction interval (solution space) L = [T^, Tf,] with the step size hx as follows 



Ti = Ta + i hx, 



0, 



Nt 



Nt 



(44) 



6.1 Line Integration Over Approximate Level Sets 

This approach is a variation of the geometric integration technique developed by Min and 
Gibou 15 , 16 . The main idea behind both methods is decomposition of the domain Q 



into simplices, which in our simplest 2D case is represented by triangulation (40), and then 
approximating the level sets given by 0(s, x) = with piecewise linear splines inside each 
simplex (triangle). While in the geometric integration approach of Min and Gibou one uses 
linear interpolation to refine locally the finite elements which contain the level sets (j){s, x) = 
and then the second-order midpoint rule for approximating line integrals over the selected 
simplices, in the present method we employ analogous approximations of the level sets, but 



without local refinement, to reduce line integral (38) to a ID definite integral which is then 
evaluated using standard quadratures. 
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The starting point for this approach is formula (38). For a fixed value of s the corre- 
sponding level set can be described as 



Mis) 

r. = U r^., (45) 

where C and M{s) is the total number of the finite elements containing segments of 
the level set F^. We thus need to approximate Lj ^ da, i.e., the line integral over 



the part of the level-set curve contained in the j-th finite element In view of (42), the 



integrand expression gis, x) = - — — -r can be approximated as 

|V0(s,x)| 

^(s,x)b^=SLif5j.^Ux), 

where q\ are the known nodal values of the function ^(s,x). An approximation F^ of the 
part of the level set F^ belonging to the j-th finite element can be obtained in an explicit 
form y = y{x), x G [x', x"], or a parametric form x = x{t), y = y{t) with t G [t', t"], based on 



representation (42) of the level-set function 0(s,x). This leads to the following two possible 



reductions of the line integral to a definite integral 

g{s,x.)da = / g{x,y{x))\ i \ +ldx (46a) 




^~g{s,^)da = 1^'^ (46b) 

which can be evaluated using standard quadratures for ID definite integrals. We then have 

fis)^j:fJ!^ fgis,^)da. (47) 

We note that the accuracy of this approach is mainly determined by the order of interpolation 
used to represent the level set F^ and the integrand expression g{s, x) which depend on the 



type of the finite elements used |33]. (The error of the quadrature employed to evaluate (46) 
does not have a dominating effect.) As was mentioned in [15], the use of triangulation (40) 
together with linear interpolation of 0(s,x) and g{s,x.) results in the overall second-order 
accuracy of this method. 

6.2 Approximation of Dirac Delta Measures 



This approach has formula (37) for its starting point and relies on a discrete approximation 



of the Dirac delta function obtained by Smereka in 11 . It is derived via truncation of 



the discrete Laplacian of the corresponding Green's function. Suppose the domain is 
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covered with a uniform Cartesian grid corresponding to (41), i.e., with nodes Xi = Xq + ih, 
Vj = yo + jh, where i,j are integer indices, xo,yo G M and h is the step size. The first-order 
accurate approximation of the discrete Dirac delta function at the node {xi,yj) is 

^(0M)=5;r^+c^^+^!r^+5;/, (48) 



where 



if < 0, 



0, otherwise, 

\(pi-i,jD%j\ 



if < 0, 



X r«J 1 1 ^ Or«J I 

0, otherwise, 



h,j+iD^(pi,j\ 



if (j)u(pij+i < 0, 



0, otherwise. 



h,j-lDy4>i,j\ 



if < 0, 

0, otherwise. 



where for the discretized level-set function 0jj = Xi, yj) we have the following definitions 



and 



— ^ ) J-^x Vi,j — ^ ) ^x'rij — 



in which e ^ 1 is used for regularization |11|. The expressions Dyipij, D~(j)ij, Dyipij are 
defined analogously. Using the definition of the discrete delta function from (48), the value 
f{s) in (37) can be thus approximated in the following way 

f{s)^h'J2kJ9^,, (49) 
where gtj are the nodal values of the function 5'(x). We note that this method was validated 



m 



11 exhibiting the theoretically predicted first order of accuracy only in cases in which 



the level sets Tg do not intersect the domain boundary d^l, a situation which may occur in 
the present reconstruction problem. 
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6.3 Approximation of Contour Integrals with Area Integrals 



Our third method, in which the level-set integral (37) is approximated with an area integral 



defined over a region containing the level set Tg, cf. (36), appears to be a new approach and 



will be presented in some detail here. It consists of the following three steps 

1. for a fixed value of the state variable s = Tj we define the interval [T^_i,T-_^i] 
[s — \hT, s + \hT\ C L; then, we have 



f{s) 



hi 



T Js-lhT 



fiOdC, 



(50) 



2. now we define a subdomain Qs,hT C Q which contains all the points of Q that lie 
between the two level-set curves T^_if^^ and F^,,,!^^ 



^s,hT — 



1, Ij. ' 

-hT,S + -hT 



X G fi, T(x) e 

see Figure [2^; we then approximate Qs,hT with the region 



(51) 



s,hT 



U where 



3 = 1 



fi? : x° G and T(x° 



see Figure|2|D, which consists of the quadrilateral finite elements f^^/^^.j, . 
with the center points x° satisfying the condition 7'(x°) G [s — (l/2)/ir, s + (l/2)/ir] 



s- ^hT,s + -Kt 
(52) 

j 1) • • • ) -^s,hT 1 



3. in view of (50), expression (37) is approximated with an area integral over the region 

which is in turn approxi- 



contained between the level-set curves F, and F 



S-^hT 



mated by the FEM region ^s,hT given by (52); finally, the integral over this region is 
approximated using the standard 2D compound midpoint rule as 



fis) 



h 



s+hhr 



T J s-\hT Jn 



9 -^s,hj-' 

n 



5(.(x)-C)^7(x)rfxrfC^^ 



(53) 



i=i 



As regards the accuracy of this approach, we have the following 



Theorem 6.1. Formula (53) is second order accurate with respect to the discretization of the 



space domain and first order accurate with respect to the discretization of the state domain, 
i.e., 

h' 



(2 ^S,hrp 



0{h') + 0{h2 



(54) 
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Proof. We start by integrating both sides of (37) over the interval ; T-_^_i] = [s — \hT] s + 
\hT\ obtaining 



s+\hT 



/(C) 



s-\hT 



5{T{^) - C) dC 



where the characteristic function 

X[T',T"](T(x)) 



(55) 



1 for T(x) G [r,T"] 
for T(x) ^ [T',T"] 



describes the subdomain ^s,hT introduced earher in ( [5l| , cf. Figure Now, using a second- 
order accurate midpoint rule for ID integration, we can express the LHS of (55) as 



■hiiT 



/(C) rfC = / ' ' ' 



■ hT + 0{hl) = f{s)hT + &{hl). 



(56) 



Approximation of the RHS in (|55|) takes place in two steps. In the first step we approximate 



the actual integration domain VLs^Ht with the union of the finite elements VLs^Kt^ cf. (52). In 
order to estimate the error 



El 



of this step, we divide the set of cells VLs^Ht iiito two subsets, see 



Figure 



where ^l^hr consists of the cells with all 4 vertices satisfying the condition T(xjt) G 

[s — ^Ht, s+^hx]. The subregion Q'^ defined as the compliment of fi* in fis^^j,, represents 
the union of "truncated" cells, i.e., cells which have at least one node outside ^s,hT- This 
subregion is in turn further subdivided into two subsets, i.e.. 



s,hT 



where 



Q' A 



s,hj^;out 



xGfi' , r(x)G 



xGli' ,T(x)^ 



-hr, s 



1. 1, 

s - -hT,s+ -hr 



We have to define one more set fi"^^ which consists of the cells with at least one vertex 
{^k}k=i satisfying the condition T(xa:) G [s — \hT,s + |/iT], but whose center points x° lie 
outside O^s^hT- We also further subdivide this set into two subsets 



O" - O" II O" 
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where 



^^;',...n= xg(]'' ,T(x)e 



1. 1. ' 

s- -hT,s+ -Kt 
s - -hT,s+ -hr 



We thus have 



5((x) dx -- 



g{x) dx + / 5f(x) dx + / (yf(x) dx, (57) 



s,n'j' 



g{x.) rfx + / g{x.) rfx 
Jn' . 

s ,nip 

g{x) (ix + / g{x) dx+ g{x) dx, 

Jn' ^ .. Jn' 



(58) 



\,hrp-,irL s,hj^ \ out 

so that the domain approximation error can be estimated as follows 



El 



fi'(x) rfx 



s.hrj^ '.out 



dx 


< 











fi'(x) (ix 



5((x) dx 



< max |5'(x) 



S,hrT 



n'' 



O" 



(59) 

where = meas Q. The second error in the approximation of the RHS of (55) is related to 
the accuracy of the quadrature applied to o(x) (ix and for the 2D compound midpoint 

rule is E2 = 0{h'^), so that we obtain 



0(/l'). 



(60) 



Comparing (53) with (56) and dividing both sides into hx we finally obtain (54) which 
completes the proof. □ 

So far, we have considered the discretizations of the physical and state spaces, f2 and L, 
as independent. We remark that using the relationship 



min |VT(x)| ■ h < Ht < max |VT(x)| ■ h 



(61) 



one could link the corresponding discretization parameters h and Ht to each other. 
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Figure 2: Illustration of approach ^^3 where line integral (37) is approximated with an area 
integral (see Section 6.3): (a) region Qs,hT which lies between the two level-set curves F^.i^^^ 

and T^^if^^ and (b) its approximation with the region flg^hr = ^Xhr ^ s,^t' "^^ere checked 

cells represent ^'^./ly shaded cells represent ^s^^y- Figure (b) also shows a part of the 

region fig^hr ~ ^"s,hT;in ^ ^'s,hT]out represented by 2 elements in the top right corner. 
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7 Computational Results 



7.1 Comparison of Different Approaches to Gradient Evaluation 

In this Section we discuss the accuracy and efficiency of the three methods for evaluation of 



expression (37) presented in Sections 6.1, 6.2 and 6.3 In order to assess their utility for the 



parameter reconstruction problem studied in this work, we will consider the following three 
test cases 

(i) single (fixed) value of s with 0(s,x) and g{x.) given analytically, 

(ii) parameter s varying over a finite range with 0(s,x) and g{x.) given analytically, 

(iii) parameter s varying over a finite range with 0(s, x) and g{x.) given in terms of solutions 
of the direct and adjoint problem. 

Tests (ii) and (iii) with s spanning the entire interval L are particularly relevant for the 
present reconstruction problem, as they help us assess the accuracy of the cost functional 
gradients over their entire domains of definition, including the values of s for which the level 
sets Ts intersect the domain boundary dQ. Results of tests (i)-(iii) are presented below. 

7.1.1 Tests for a Single Value of s with 0(s,x) and g(x) Given Analytically 

Here we employ our three methods to compute numerically the value of a line integral over 
the circle x"^ + y'^ = 1 

(a) in domain Qi = [—2, 2]^ which contains the entire curve 

Ie.,1 = [ (3x2 _ y2^ ^ 27T (62) 



(this test problem is actually borrowed from [TT]), 
(b) and in domain Q2 = [0, 2]^ which contains only a part of the curve in the ffist quadrant 



Iex,2= I (3x2-y2)rf^=^. (63) 



The main difference between test cases (a) and (b) is that while in (a) the contour is entirely 
contained in the domain Qi, it intersects the domain boundary 8^2 in case (b). As shown 
in Figure [3| methods #1 (line integration) and ^^3 (area integration) in both cases show 
the expected accuracy of 0{h'^), where h = Ax = Ay = 2~'^^+*\ i = 1 ... 6, while method 
#2 (delta function approximation) is 0{h^^'^) accurate in case (a) and only 0(h^) accurate 
in case (b). We also add that the line integration method exhibits the smallest constant 
prefactor characterizing the error. 
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(a) 



(b) 



Figure 3: Relative error |-^ — 1\, i = 1,2, versus discretization step h = Ax = Ay in the 
numerical approximation Ih of (a) line integral (62) and (b) line integral (63), see Section 



7.1.1 Triangles represent the line integration approach (method #1), circles represent the 
results obtained using the delta function approximation (method #2), whereas asterisks 
show the data for the area integration approach (method #3). 



7.1.2 Tests for s Varying Over a Finite Range with (j){s,x) and g{x) Given Ana- 
lytically 

In order to analyze this case we will introduce a new diagnostic quantity. We begin with the 



integral transform formula (28) applied to some perturbation yu'(T(x)) 



/i'(r(x)) 



+ 00 



5(0(s, x))yu'(s) ds, 



(64) 



where 0(s,x) = T(x) — s. Multiplying both sides of (64) by 5'(x), integrating over the 



domain fl and changing the order of integration we obtain the following useful identity 



/w'(^(x))fl'(x) rfx 



+00 



/(s)/i'(s) ds 



(65) 



with /(s) defined in (37), where the RHS has the structure of the Riesz identity for the 



Gateaux differential of the cost functional, cf. (30), whereas the LHS is a simple area integral 
which can be easily evaluated using high-accuracy quadratures. Given the formal similarity 



of the RHS of (65) and the Riesz formula (30), this test is quite relevant for the optimization 



problem we are interested in here. We will thus use our three methods to evaluate the RHS 



of (65) and compared it to the LHS, which is evaluated with high precision on a refined grid 



in n, so that it can be considered "exact". Our tests are based on the following data 



24 



• spatial domain f2 = [0, 1]^ discretized with h = Ax = Ay = 2 (^+*), i = 1 . . . 7, 

• state domain L = [Ta,Tb], where Ta = 100, Tf, = 700, discretized using Nt = 
200, 1000, 10000 points for methods #2 and #3 and Nt = 200, 2000, 20000 points 

for method #1: hx = 

Nt 

• T(x) = 100(a;2 + ^ 390, ^(x) = cos(x) + 3sin(2?/ - 1), x e Q, 

• perturbations used n[{s) = exp(-j^), n'^is) = ^ and n'^i-s) = + ^ + h 

As is evident from Figure |4| all three methods show similar qualitative behavior, namely, 
the error decreases with decreasing h until it saturates which is due to the error terms de- 
pending on Ht becoming dominant. The saturation value of the error depends on the state 
space resolution and is different for the different methods. Method 7^8 (area integration) 
reveals accuracy 0{h'^), whereas method 7^2 (delta function approximation) is again only of 
accuracy about 0{h) for the same discretization of the interval L. Method #1 (line integra- 
tion) performs better and shows accuracy up to 0{h^), but requires much finer resolution in 
the state space, namely Nt > 20000 (/it < 0.03) is needed for this behavior to be visible. On 
the other hand, method #3 (area integration) leads to the smallest errors for all the cases 
tested. 

Analogous data is plotted in Figure [5] now as a function of the state space resolution Ht 
with h acting as a parameter. Similar trends are observed as in Figure |4| namely, the errors 
decrease with Ht until they eventually saturate when the error terms depending on h become 
dominant. Methods #1 and #2 reveal accuracy 0{hT), whereas method #3 has accuracy 



0(/iy^~^) which is actually better than stipulated by Theorem 6.1, cf. (54). Method #3 is 
also characterized by the smallest value of the constant prefactor leading to the smallest 
overall errors. 



7.1.3 



Tests for s Varying Over a Finite Range with 
Solutions of Direct and Adjoint Problem 



(s,x) and g(x) Given by 



s,x) = T(x) — s and g{x.) = 
where the fields u, T and u* come from the solutions of 



usmg 



We now repeat the test described in Section 7.1.2 
[Vu(x) + (Vu(x))^] : Vu*(: 
the direct and adjoint problem ([T|-(|2]) and (18)-(19) at some fixed time t, see Figure |6 
(details of these computations will be given in Section 7.5). As before, we discretize the 



domain = [0, 1]^ with the step h = Ax = Ay = 2 i = 1 ... 7, and the state space 

L = [Ta,Tb\, Ta = 100 and = 700, with the step Ht = 0.06 {Nt = 10000). 

The data shown in Figure [7] confirms our findings from Sections |7. 1. 1| and 7.1.2 , namely, 
that in this case as well the error of all three methods decreases with h until it eventually 
saturates when the errors depending on Ht become dominant. Method #3 is again char- 
acterized by the smallest prefactor and hence leads to much smaller overall errors than in 
methods 7^1 and 7^2. The computational complexity of our three approaches is addressed 
in Figure [8| where Ne is defined as a number of computational elements, i.e., Ne = Na, or 
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Figure 4: Relative error |^ — 1|, where ^ 
in approximating the RHS in (65). 



RHS of lesi "' '^^^sus discretization step h = Ax = Ay 
The first, second and third row of figures show the results 
for ji^ and /^g, respectively, while the figures in the first, second and third column represent 
line integration (#1), delta function approximation (#2) and area integration (#3) methods, 
respectively. 
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(g) 



(h) 



- N = 1 28 
-N = 512 

- N = 2048 




(i) 



Figure 5: Relative error |^ — 1|, where ^ 
state space L in approximating the RHS in (|65|). 



RHS of lesi ' '^srsus discretization step Ht in the 
e first, second and third row of figures 



show the results for /x'^, ^2 and n'^, respectively, while the figures in the first, second and 
third column represent line integration (#1), delta function approximation (#2) and area 
integration (#3) methods, respectively. 
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Ne = Nu, using (40) or (41), respectively. We see that, while the complexity of methods #1 
and #2 scales as @{^/N~e), method #3 exhibits the scaling of 0(A^e)- On the other hand, 
however, method ^^3 has the smallest prefactor and, at least in the range of resolutions 
considered here, results in the shortest execution time. 

In conclusion, these observations make the area integration approach (method #3) the 
method of choice for the present parameter reconstruction problem, and this is the approach 
we will use in all subsequent computations. 



7.2 Models for Constitutive Relations 

For validation purposes one needs an algebraic expression to represent the dependence of 
the viscosity coefficient on temperature which could serve as the "true" material property 
we will seek to reconstruct. The dynamic viscosity in liquids is usually approximated by 
exponential relations |34| and one of the most common expression for the coefficient of the 
dynamic viscosity is the law of Andrade (also referred to as the Nahme law) which is given 



in the dimensional form valid for T expressed in Kelvins in ( 66 ) below 

fi{T) = Cie^^/^, (66) 

where Ci, C2 > are constant parameters. As regards the thermal conductivity /c, since it 
typically reveals a rather weak dependence on the temperature, for the sake of simplicity we 
will treat it as a constant setting k = 0.002 in all computations presented in this paper. 
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Figure 7: Relative error |^ — 1|, where ^ 



LHS of J65l 



Ax 



= Ay 
Vu* 



RHS of (l65li ' ^^^sus discretization step h 
in estimating the RHS in ( [65|), where 0(s,x) = T(x) — s, and g{x) = [Vu + (Vu)-^] 
is obtained by solving ([T|-(|2| and (18)-(19). Figures (a), (b) and (c) show the results 
for ^2 ci-iid /i3, respectively, using the same discretization of the state space L with 
Nt = 10000. Triangles represent the line integration approach (#1), circles show the results 
for the method of the delta function approximation (#2), while asterisks show the data from 
the area integration approach (#3). 




Figure 8: CPU time (in seconds) versus the number Ne of computational elements used in the 
different approaches, namely, (triangles) finite elements f2f in the line integration approach 
(method #1), (circles) grid nodes in the method of the delta function approximation (#2) 
and (asterisks) finite elements in the area integration approach (#3). The data shown 
corresponds to the estimation of the RHS in (65) for ^[ using the same discretization of the 
state space with Nt = 10000, 0(s, x) = T(x) — s and (^(x) = [Vu + (Vu)-^] : Vu* obtained 
by solving and (|l8|)-(|l9|). 
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Figure 9: Geometry of the 2D lid-driven (shear-driven) cavity. 
7.3 Model Geometry and PDE Solvers 

To validate the accuracy and performance of the proposed approach to reconstruct /i(T), we 
use a simple 2D lid-driven (shear-driven) cavity flow, cf. Figure [oj as our model problem. 
Due to the simplicity of its geometry and boundary conditions, the lid-driven cavity flow 



problem has been used for a long time to validate novel solution approaches and codes 35 36 
Numerical results are available for different aspect ratios and the problem was solved in both 
laminar and turbulent regimes using different numerical techniques. Thus, this problem is a 
useful testbed as there is a great deal of numerical data that can be used for comparison. Our 
code for solving direct problem ([T])-([2]) and adjoint problem (18)-(19) has been implemented 



using FreeFem++ [37j, an open-source, high-level integrated development environment for 
the numerical solution of PDEs based on the the Finite Element Method. The direct solver 
has been thoroughly validated against available benchmark data from |35l|36| , and all details 



are reported in 31 



To solve numerically the direct problem we discretize system ([T])-([2]) in time using a 
second-order accurate semi-implicit approach. Spatial discretization is carried out using 



triangular finite elements (40) and the P2 piecewise quadratic (continuous) representations 
for the velocity u and the temperature T fields, and the PI piecewise linear (continuous) 
representation for the pressure p field. The system of algebraic equations obtained after such 
discretization is solved at every time step with UMFPACK, a solver for nonsymmetric sparse 
linear systems |38|. We add that incompressibility is ensured by an implicit treatment of 



equation (lb). Stability is enforced by choosing the time step At so that it satisfies the 
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following CFL condition 



At 



< 1. 



(67) 



The same technique is used for the numerical solution of adjoint problem (18)-(19). 

All our computations are performed using a 2D square domain Q = [0, 1]^ shown in Fig- 
ure |9} Governing system ([T])-([2]) and adjoint system (18)-(19) are discretized on a uniform 



mesh with N = N^, = Ny = 32 grid points in every direction using triangular finite elements 
combined with the cubic spline interpolation of the function /x(T(x)). The rather modest 
spatial resolution used is a consequence of the fact that in a single reconstruction problem 
the governing and adjoint systems need to be solved 0(10^ — 10"^) times, hence there are lim- 
itations of the computational time. Unless stated otherwise, the interval L = [100.0, 700.0] is 
discretized using an equispaced grid with = 600 points. The actual constitutive relation 
/i(T) we seek to reconstruct is given by Andrade law (66) with Ci = 0.001 and C2 = 10^. In 
the computational tests reported below we used M = 9 measurement points distributed uni- 
formly inside the cavity (Figure [9]). To mimic an actual experimental procedure, first relation 
(66) is used in combination with governing system ([T])-([2]) to obtain pointwise temperature 
measurements {Tj}^]^. Relation (66) is then "forgotten" and is reconstructed using gradient- 



based algorithm ([14j). In terms of the initial guess in (14), unless stated otherwise, we take a 
constant approximation /zq to (66), given by /iq = ^ (/^(^a) + /^(^/s)) = ^ (^e'^^^^°' + e'- 
which translates into the following expression for the new optimization variable 6, cf. 
6*0 = a/z^o — "^M' where = ^^fxiTp) = ^e'^'^^'^^. Since in the present problem the viscosity 
/x(T) is a function of the temperature, the Reynolds number is defined locally (both in space 
and in time) and varies in the range Re = 0.05 -^ 240. 



Unless stated otherwise, the boundary conditions for the temperature are Tg 



top 



500 and Tb 



else 



Ta = 300 which results in the identifiability region I = [300.0,500.0]. 



The velocity boundary conditions = [ub,vb] are given by wsltop = ^Tq cos(27rt), Uq = I 
and ffiltop = on the top boundary segment and u^leise = on the remaining boundary 
segments. Their time-dependent character ensures that the obtained flow is unsteady at the 
values of the Reynolds number for which self-sustained oscillations do not spontaneously 
occur (the study of higher Reynolds numbers was restricted by the numerical resolution used, 
see comments above). The initial conditions {uo,To} used in the reconstruction problem 
correspond to a developed flow obtained at t = 10 from the following initial and boundary 
conditions T^ltop = 500, T^leise = 300 and usltop = 1, ^^sltop = 0, u^leise = 0, Tq = 300, 
Uq = 0. We emphasize that adjoint system (18)-(19) is in fact a terminal-value problem 



which needs to be integrated backwards in time, and its coefficients depend on the solution 
{u, T} of the direct problem around which linearization is performed at the given iteration. 
Our reconstructions are performed using the following time windows [0,t/], tf = ^, l} 
which correspond to a fraction of, or a full, forcing cycle in the boundary conditions described 
above. These time windows are all discretized with the time step At = 5 ■ 10~^ in both the 
direct and adjoint problems. This choice of the time step At ensures stability by satisfying 



the CFL condition (67) 
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7.4 Validation of Gradients 



In this Section we present results demonstrating the consistency of the cost functional gra- 
dients obtained with the approach described in Section |4} In Figure [lO] we present the L2 
and several Sobolev gradients obtained at the first iteration. In the first place, we ob- 
serve that, as was anticipated in Section |4], the L2 gradients indeed exhibit quite irregular 
behaviour lacking necessary smoothness which makes them unsuitable for the reconstruction 



of constitutive relations with required properties, cf. (10). On the other hand, the gradients 



extracted in the Sobolev space are characterized by the required smoothness and there- 



fore hereafter we will solely use the Sobolev gradients. Next, in Figure 11^ we present the 
results of a diagnostic test commonly employed to verify the correctness of the cost func- 
tional gradient It consists in computing the directional Gateaux differential J'{0]0') 
for some selected perturbations 6' in two different ways, namely, using a finite-difference 



approximation and using (29) which is based on the adjoint field, and then examining the 



ratio of the two quantities, i.e.. 



,-1 



[j{e + tO') - j{d)] 

j^^VeJ{s)9'{s)ds 



(68) 



for a range of values of e. If the gradient ^eJ^iO) is computed correctly, then for intermediate 
values of e, K{e) will be close to the unity. Remarkably, this behavior can be observed 



in Figure 11 over a range of e spanning about 6 orders of magnitude for three different 
perturbations 6'{T). Furthermore, we also emphasize that refining the time step At used 
in the time-discretization of ([T])-([2]) and (18)-(19) yields values of ^(e) closer to the unity. 
The reason is that in the "optimize-then-discretize" paradigm adopted here such refinement 



of discretization leads to a better approximation of the continuous gradient (29). We add 



that the quantity log;^o l'^(^) ~ M plotted in Figure 11 d shows how many significant digits of 
accuracy are captured in a given gradient evaluation. As can be expected, the quantity ^(e) 
deviates from the unity for very small values of e, which is due to the subtractive cancellation 
(round-off) errors, and also for large values of e, which is due to the truncation errors, both 
of which are well-known effects. 



7.5 Reconstruction Results 



We solve minimization problem (|9| using the Steepest Descent (SD), Conjugate Gradient 

(CG) and BFGS algorithms |26| and, unless indicated otherwise, using Sobolev gradients 

computed with i = 200.0 which was found by trial-and-error to result in the fastest rate of 

j'(e("))-j'(6>("-^)) 



convergence of iterations (14). The termination condition used was 
The behavior of the cost functional J'{9^^^ 



< 10 



-6 



as a function of the iteration count n is shown in 
Figure 121 for all three minimization algorithms (SD, CG and BFGS). We note that in all 
cases a decrease over several orders of magnitude is observed in just a few iterations. Since the 
three descent methods tested reveal comparable performance, our subsequent computations 
will be based on the steepest descent method as the simplest one of the three. The effect 
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Figure 10: Comparison of (thin solid line) the L2 gradient 'Vg^J' and the Sobolev gradients 
defined in (31) for different values of the smoothing coefficient (thick dashed line) 
i = 2.5, (thick dash-dotted line) i = 10.0 and (thick sohd line) i = 200.0 at the first 
iteration with the initial guess /xq = const = 0.0177. The vertical dashed lines represent the 
boundaries of the identifiability interval I and the vertical scale in the plot is arbitrary. 
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Figure 11: The behavior of (a) ^(e) and (b) log^^g l'^(^) — 1| as a function of e for different 
perturbations (triangles) 9'(T) = ^, (circles) 9'{T) = e^wm and (asterisks) 9'{T) = —(^^ + 
^ + |. The time steps used in the time integration of (jT])-(|2| and (18)-(19) are (dash-dotted 
line) At = 5.0 ■ 10"^ and (solid line) At = 5.0 • 10"^ 



of the different initial guesses /iq on the decrease of the cost functional is illustrated in 
Figure [T2]d. Again, a substantial decrease of the cost functional corresponding to about 5-6 
orders of magnitude is observed for all the initial guesses tested. Reconstructions fi(T) of 
the constitutive relation obtained using the initial guess /xo = ^ (^^C2/Ta _j_ gC2/T0j _ 0.0177 
and the optimization time windows with tf = are shown in Figure 13 Comparing the 



1 1 

4' 2' 

accuracy of the reconstruction obtained for these different time windows, we can conclude 
that better results are achieved on shorter time windows tf = |, ^. Given considerations of 
the computational time, hereafter we will therefore focus on the case with tf = l. In Fi gure 



14 we show the reconstructions fi{T) of the constitutive relation obtained from different 



initial guesses already tested in Figure 12 d such as constant values of /xq, I^o{T) varying 
linearly with the temperature T and /iq given as a rescaling of the true relationship /t. As 
expected, the best results are obtained in the cases where some prior information about the 
true material property is already encoded in the initial guess /iq, such as the slope, cf. Figure 
14:;, or the exponent, cf. Figure [T4|i. We may also conclude that, since all the reconstructions 



shown in Figures [13| and [14| are rather different, the iterations starting from different initial 
guesses converge in fact to different local minimizers. However, it should be emphasized 
that in all cases the reconstructions do capture the main qualitative features of the actual 
material property with differences becoming apparent only outside the measurement span 
interval Ai. In order to make our tests more challenging, in the subsequent computations 



13) which contains 



we will use the initial guess /xq = | (/i(^a) + MTp)) = 0.0177 (cf. Figure 
little prior information about the true material property. 

In Figure 15 we show the results of the reconstruction performed with a larger identifia- 
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Figure 12: (a) Decrease of the cost functional J'{9^'^^) with iterations n using the Sobolev 
gradient defined in (31) and obtained with (sohd hne) steepest descent, (dash dotted 

hne) conjugate gradient and (hne with dots) BFGS methods with initial guess /io = 0.0177. 
(b) Decrease of the cost functional J {6^'^^) with iterations n for different initial guesses: 
(dots) /io = p^iTa) = 0.0280, (dash-dotted line) /iq = /i(71) = 0.0042, (dashed line) fio(T) 
varying linearly between jl{Ta) and jl{Tp), (thin solid line) /io ~ ^ ~' 
/io = 0.0177. 



/i.(T) and (thick solid line) 



bility region X = [250.0, 700.0] = L which is done by adopting the corresponding values for 
temperature boundary conditions (2b), i.e., Tb |eise = 250 and Tgltop = 700. We note that in 



this case the target interval for the reconstruction L has been chosen to coincide with iden- 
tifiability interval X. Fairly accurate reconstruction can be observed in this problem as well, 
and we emphasize that this is also the case for values of the temperature outside the iden- 
tifiability interval studied in the previous case. This example demonstrates that accurate 
reconstruction on different intervals L can in fact be achieved by adjusting the identifia- 



bility region via a suitable choice of temperature boundary conditions (2b). This process 
can be interpreted as adjusting the conditions of the actual experiment used to obtain the 
measurements {Ti}^^. 

7.6 Reconstruction in the Presence of Noise 

In this Section we first assess the effect of noise on the reconstruction without the Tikhonov 
regularization and then study the efficiency of the regularization techniques introduced in 
Section [5] In Figure [l6^,b we revisit the case already presented in Figure [l3^ (reconstruction 
on the interval L = [100.0,700.0] with the identifiability region X = [300.0,500.0]), now for 
measurements contaminated with 0.05%, 0.1%, 0.3%, 0.5% and 1.0% uniformly distributed 
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Figure 13: Reconstruction fi{T) of the material property obtained using the Sobolev gra- 



dients defined in (31) on (a,b,c) the interval L and (d) close-up view showing the interval 
outside the identifiability region X with the time window [0,t/], where (a) tf = j (b) t/ = ^ 



and (c) tf = 1. The dash-dotted line represents the true material property (66), the thick 
solid, dashed and dash dotted lines are the reconstructions for (a,d) = |, (b,d) t/ = | and 
(c,d) tf = 1, respectively, whereas the dashed line represents the initial guess /io = 0.0177; 
the vertical dash-dotted and dotted lines represent, respectively, the boundaries of the iden- 
tifiability interval X and the measurement span A4. 
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Figure 14: Reconstruction fi{T) of the material property obtained using different initial 
guesses (a) /xq = P^iTa) = 0.0280, (b) /iq = f^iTb) = 0.0042, (c) fio{T) varying linearly 



between n{Ta) and Ai(T^) and (d) fio = ^fi{T), and the Sobolev gradients defined in (31) 



on the interval L. The dash-dotted line represents the true material property (66), the 
solid line is the reconstruction fi{T), whereas the dashed line represents the initial guess 
/io; the vertical dash-dotted and dotted lines represent, respectively, the boundaries of the 
identifiability interval X and the measurement span M.. 
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Figure 15: Reconstruction fi(T) of the material property using an extended identifiability 
region X = [250.0, 700.0] = L shown on (a) the interval [100.0, 700.0] and (b) magnification 
of this new identifiability region. The dash-dotted line represents the true material property 
(66), the solid line is the reconstruction fi{T), whereas the dashed line represents the initial 
guess /io; the vertical dash-dotted and dotted lines represent, respectively, the boundaries 
of the identifiability interval Xq = [300.0, 500.0] used previously and the measurement span 
M, while the dashed vertical lines show the boundaries of the new identifiability interval X. 
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noise and without Tikhonov regularization. To incorporate noise, say of •)]%, into the mea- 
surements {Tj(t)}*f]^, we replace these measurements at every discrete time step tj G [0,t/] 
with a new set {Ti{tj)}fii, where the independent random variables T^{tj) have a uniform 
distribution with the mean Ti{tj) and the standard deviation Ar] = -^Yl!iLi'^i{^j) ' Too%- 
Unless stated otherwise, in order to be able to directly compare reconstructions from noisy 
measurements with different noise levels, the same noise realization is used after rescaling to 
the standard deviation Ar]. As expected, in Figure 16i,b we see that increasing the level of 
noise leads to oscillatory instabilities developing in the reconstructed constitutive relations 
fi{T). We note that the reconstructions become poor already for relatively low noise levels, 
i.e., on the order of 1%. One reason for this seems to be the time-dependent nature of the 
problem in which independent noise is added to the measurements at every (discrete) time 
instant leading to accumulation of uncertainty. Indeed, such loss of information was not 
observed in the case of the steady problem studied in [7j where reliable reconstructions could 
be obtained with noise levels an order of magnitude larger than in the present investigation. 
In regard to the results shown in Figure [T6|i,b, we add that the pattern introduced by the 
noise in the reconstructions depends on the specific noise sample used (which was the same 
for all the reconstructions shown in the Figure). Reconstructions performed using different 
realizations of the noise produce distinct patterns in the reconstructed constitutive relations 
fi{T). In our computational experiments we also observe that inclusion of noise in the mea- 
surements tends to replace the original minimizers with perturbed ones (this is evident from 
the uniform, with respect to T, convergence of the perturbed minimizers to the noise-free 
reconstructions as the noise level goes to zero in Figure 16i,b). 

The effect of the Tikhonov regularization is studied in Figure Il6b,d, where we illustrate 



the performance of the technique described in Section |5| cf. (33), on the reconstruction 
problem with 1.0% noise in the measurement data (i.e., the "extreme" case presented in 
Figures 16i,b). In terms of the (constant) reference function we take 9 = y/fio — ^fi, where 
IIq = 0.0177. We note that by increasing the values of the regularization parameter A in 
(33) from (no regularization) to 2500 we manage to eliminate the instabilities caused by 
the presence of noise in the measurements and obtain as a result smoother constitutive 



relations, cf. Figure [T6p,d. We add that, while after introducing the Tikhonov regularization 
the reconstructed solutions converge in fact to different local minimizers (in comparison 
with the reconstructions without noise), this does not prevent the reconstructions from 
capturing the main qualitative features of the actual material property. Systematic methods 
for determining the optimal values of regularization parameters are discussed for instance in 



30 . Finally, in Figure 17 we present the relative reconstruction errors — /i||Li(x) / ||/i||Li{x) 
obtained using the approach described earlier in Section [5] for data with different noise levels 
and averaged over 10 different noise samples. From Figure [TT] we conclude that larger values 
of the regularization parameter A are required for more noisy measurements. We close this 
Section by concluding, in agreement with our earlier results reported in |7], that Tikhonov 
regularization performs as expected in problems with noise present in the measurement data. 
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Figure 16: (a,b) Reconstruction fi{T) of the material property obtained in the presence 
of different noise levels in the measurement data: (thick solid line) no noise, (dotted line) 
0.05%, (dashed line) 0.1%, (dash-dotted line) 0.3%, (thin solid line) 0.5%, and (thick dashed 
line) 1.0% on (a) the interval L and (b) close-up view showing the identifiability interval 
X. (c,d) Effect of Tikhonov regularization on the reconstruction from the measurement 
data with 1.0% noise using regularization term (33) on (c) the interval L and (d) close-up 
view showing the identifiability interval X. In both figures (c,d) the following values of the 
regularization parameter were used: (thick dashed line) A = 0, (circles) A = 2.5, (dashed line) 
A = 6.25, (thin solid line) A = 25.0, (dash-dotted line) A = 250.0, and (dots) A = 2500.0. 
For all figures the horizontal dashed line represents the initial guess /io = 0.0177; the vertical 
dash-dotted and dotted lines represent, respectively, the boundaries of the identifiability 
interval X and the measurement span Ai. 
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Figure 17: Relative Li reconstruction errors — /i||Li(x) / ||/i||Li(x) obtained in the presence 
of noise with the amphtude indicated and averaged over 10 samples: (dash-dotted line) 
reconstruction with Sobolev gradients and without Tikhonov regularization, and (solid line) 
reconstruction with Tikhonov regularization term ([33)) [(circles) A = 2.5, (triangles) 



A = 250.0]. The thick dashed line represents the "error" in the exact material property (66) 
obtained by adding noise to T and averaging over time steps. 



8 Conclusions and Summary 

We developed an optimization-based approach to the problem of reconstructing temperature- 
dependent material properties in complex thermo-fluid systems. As a model problem we 
considered two-dimensional unsteady flows in a lid-driven cavity involving also heat trans- 
fer. A key element of this approach is the gradient of the cost functional which is computed 
based on a suitably-defined adjoint system, and is shown to have mathematical structure 
different than in most problems involving adjoint-based PDE optimization. We discussed 
three different numerical approaches to evaluation of these cost functional gradients which 
are given in terms of integrals defined on the level sets of the temperature field. As compared 
to earlier work on the numerical approximation of such expressions , we also addressed 



at length the question of the discretization of the solution space which is specific to our 
reconstruction problem. Evidence is shown for the superior performance with respect to the 
discretizations of the physical and the solution space, as well as the computational time, of 
a new approach to evaluation of gradients which is proposed in this study. 

The reconstruction results obtained demonstrate good performance of the algorithm, in 
particular, it is shown that by suitably adjusting the boundary conditions in the governing 
system we can extend the identifiability region. There are two comments we wish to make 
concerning these results. As regards the data shown in Figure [T3| it might appear somewhat 
surprising that reconstructions on longer time windows (which use more information from 
the measurements) do not necessarily lead to better results. A probable reason is that 
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such reconstruction problems defined on longer time windows tend to be harder to solve 
numerically, hence the solutions found may not be global minimizers (or even "good" local 
ones). Another comment we have concerns the relatively modest noise levels for which 
stable reconstructions could be obtained in Section 17.61 We note that inclusion of even 
a small amount of noise could alter the nature of the reconstructed constitutive relations. 
On the other hand, we also remark that reconstructions performed based on a steady-state 
problem and reported in |7] did allow for significantly higher noise levels. We therefore 
conjecture that the sensitivity to noise in the present problem is related to its time-dependent 
character, as here the effect of instantaneously small noise levels may be compounded by its 
continuous accumulation over the entire time window. Since the flow problem studied here 
was admittedly very simple, this issue may certainly limit the applicability of the proposed 
method to problems characterized by higher Reynolds numbers and therefore merits a more 
thorough study in the future (this is not the question of the numerical resolution alone, 
but rather of the interplay between this resolution and the complexity of the underlying 
optimization problem). Experience with other data assimilation problems may suggest that 
acquiring the measurements less frequently in time may actually help mitigate this effect. 

The approach developed in the present study was recently successfully applied to the 
problem of system identification involving a dynamical system in the phase-space represen- 
tation. More precisely, in ^O] we used this method to reconstruct an invariant manifold 
in a realistic reduced-order model of a hydrodynamic instability. Our future work will in- 
volve extensions of the present approach to more complicated problems involving systems 
of coupled PDEs depending on time and defined on domains in three dimensions. In the 
context of such systems an interesting issue is the reconstruction of anisotropic constitutive 
relations. A more challenging problem is related to reconstruction of constitutive relations 
in systems involving phase changes. In addition to the governing PDE in the form of a 
free-boundary problem, one would also have to deal with constitutive relations with distinct 
functional forms in each of the phases. Interesting questions also arise in the reconstruc- 
tion of material properties defined on the interfaces between different phases, such as for 
example the temperature-dependent surface tension coefficient playing an important role in 
Marangoni flows. In the study of such more complicated problems close attention will need 
to be paid to the question of ensuring consistency of the reconstructed constitutive relation 
with the second principle of thermodynamics, which can be done by including a form of 



the Clausius-Duhem inequality 20 in the optimization formulation. In such more general 
problems it is not obvious whether this additional constraint can be eliminated by introduc- 
ing a slack-variable formulation similar to the one used in this study, and one may have to 
perform numerical optimization in the presence of inequality constraints. These questions 
are left for future research. 
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